Electronic thermal transport in strongly correlated multilayered nanostructures 
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The formalism for a linear-response many-body treatment of the electronic contributions to ther- 
mal transport is developed for multilayered nanostructures. By properly determining the local 
heat-current operator, it is possible to show that the Jonson-Mahan theorem for the bulk can be 
extended to inhomogeneous problems, so the various thermal-transport coefficient integrands are 
related by powers of frequency (including all effects of vertex corrections when appropriate). We 
illustrate how to use this formalism by showing how it applies to measurements of the Peltier effect, 
the Seebeck effect, and the thermal conductance. 
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I. INTRODUCTION 

As materials and device growth techniques mature and 
improve, it becomes more possible to create artificial sys- 
tems, composed of well-defined numbers of flat planes of 
one material grown on top of another material. A de- 
vice can be engineered by determining the different kinds 
of materials to grow and the thicknesses of the differ- 
ent multilayers. If we assume the growth process is per- 
fect, so the planes are atomically flat, with no interface 
roughness, then we have an inhomogeneous quantum- 
mechanical problem to solve for the behavior of electrons 
in the system, where the inhomogeneity lies in one di- 
mension only. 

The properties of devices grown in this fashion are 
more complicated if some or all of the materials that 
make up the device are composed of strongly correlated 
electron materials, where the properties of the electrons 
cannot be described solely by an independent-particle 
picture like band theory. These systems are increasing 
in interest because bulk strongly correlated materials ex- 
hibit exotic phenomena, and show promise in demon- 
strating high tunability of their properties. What is less 
studied is how these properties can be modified by con- 
finement/deconfinement effects that are possible in mul- 
tilayered nanostructures. 

In addition, there has been little theoretical develop- 
ment of the thermal transport effects in such multilayered 
systems. Some evidence from examining the interface be- 
tween two materials, indicates that thermal transport ef- 
fects in inhomogeneous systems can create large enhance- 
ments to the performancei, but the theory has not been 
fully developed within a Kubo-like context which allows 
the many-body aspects to be treated fully. Semiclassical 
approaches have also been employed^, but that theory 
is also inadequate to treat strongly correlated materials. 



Finally, the important problem of phonon transport in 
such systems has been examined extensively, but is be- 
yond the scope of what we cover in this work. 

One of the most important results in bulk thermal 
transport is the Jonson-Mahan theorem2i4, which pro- 
vides an exact relationship between different thermal 
transport coefficients. Using the Jonson-Mahan theorem 
makes calculation of the thermal transport only slightly 
more complicated than the calculation of the charge 
transport, and has enabled much of the theoretical work 
in strongly correlated thermoelectrics. Here we show how 
to generalize the Jonson-Mahan theorem to multilayered 
nanostructures, which also greatly simplifies the calcula- 
tion of the thermal transport coefficients. 

The idea to use multilayered nanostructures, or more 
complicated geometries, for enhancing thermoelectric 
performance of refrigerators or power generators was 
proposed^ in the 1990s and enhancements have been seen 
recentlj»&. The focus in that work was along the ideas 
of an electron-crystal-phonon-glass approach, where the 
nanostructures are engineered to preserve the electronic 
properties, while making the phonon transport similar 
to that in a disordered glass. It is possible that one can 
actually employ the nanostructure engineering to pro- 
duce enhancements to the electronic transport proper- 
ties, while simultaneously reducing the phonon thermal 
transport, so this basic approach may be pushed further 
than theorized in the original presentations. One key to 
being able to enhance the electronic properties, is to be 
able to tune the electron correlation properties with a 
proper engineering of the nanostructure and to engineer 
the charge redistribution at the interfaces. Before initi- 
ating such a program, one needs to be able to properly 
calculate the thermal transport in a strongly correlated 
device, and we derive the formalism for how to do this 
here. 

The systems we describe in this work involve multilay- 
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ered devices constructed of atomically flat planes which 
can be composed of different materials. Each system is 
inhomogeneous in the z-direction, which is the direction 
where the planes are stacked. We take the left lead to 
be identical to the right lead, so the system will have a 
"mirror symmetry" , and the device will have its chemical 
potential determined by that of the bulk leads. We use 
Roman letters {i, j, ...) to denote the lattice sites within 
each plane {i. e., the x — y coordinates), and Greek let- 
ters (a, /3, ...) to denote the individual planes (i.e., the 
z-coordinate). We require the system to be translation- 
ally invariant within each plane, and for the lattice struc- 
ture of each plane to be identical, so that the connection 
between planes is between corresponding sites in the two 
planes, and is the same for each site. The latter require- 
ment is by no means necessary, but it greatly simplifies 
the notation for the formalism, so we adopt it here. 

We will consider three different types of Hamiltonians 
here: (i) the Hubbard model?; (ii) the Falicov-Kimball 
model^; and (iii) the periodic Anderson model^. We use 
a multiple index ai to denote the ith planar site on plane 
a. In the Hubbard model, we have conduction electrons, 
whose creation and annihilation operators are denoted 
cj,j^ and c^jcTJ respectively, for electrons sitting at the 
lattice site denoted by ai and with z-component of spin 
a. The Falicov-Kimball model has two kinds of elec- 
trons: conduction electrons (which are described by sim- 
ilar operators as in the Hubbard model, but without spin, 
for simplicity) and localized electrons (also chosen to be 
spinless and created or destroyed by the operators /^^ 
and /Q,j). The periodic Anderson model has spin-one- 
half conduction and /-electrons, which are denoted by 
the familiar operators, except now all operators will have 
spin labels. All models can be expressed as the sum of 
two terms in the Hamiltonian — an inhomogeneous hop- 
ping term and an interaction-hybridization term. The 
inhomogeneous hopping term is essentially the same for 
all three models. It is 

Whop = — ^aii'^ia-^ajcr 
a zj'G plane o" 

^ tga+l'^'Lia'^a+litT 
a zGplanc cr 

^ ^QCt+l'^i+liCT'^Qicr (1) 

a zGplanc cr 

where we do not include a sum over spin for the spinless 
Falicov-Kimball model. We assume the hopping matri- 
ces are real symmetric matrices, and one should note that 
the hopping between planes is only between neighboring 
planes and between corresponding sites within the two 
planes. The magnitudes of the hopping matrices within 
the planes and between the planes can vary, but the pla- 
nar hopping matrices must be translationally invariant to 
go to a mixed momentum-space-real-space basis, which 
is commonly done in these types of problems. If the 
planes are square-lattice planes, then the underlying lat- 
tice topology will be that of a simple cubic lattice, but 



the hopping need not be the same everywhere. 

The interaction-hybridization term is different for each 
model. For the Hubbard model it is 

i(Eplanc 

for the spinless Falicov-Kimball model it is 
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and for the periodic Anderson model it is 

a iGplaiic cr 

+ ^ 5Z ^O'-fail fail fail fail 
a plane 
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For the Falicov-Kimball model, we often replace the term 
faifai by the symbol Wai which equals if no localized 
electrons are at site ai and equals 1 if a localized electron 
is at site ai. All interaction and hybridization terms can 
vary from plane to plane, but they must be the same for 
every lattice site within the planes to preserve transla- 
tional invariance within the planes. The total Hamilto- 
nian is then 



'H — Hhop + 'Hint - M-^i 



(5) 



for all of the models. The symbol /i is the chemical po- 
tential, and A/" denotes the electron number operator, 
chosen to be the total conduction electron number op- 
erator for the Hubbard and Falicov-Kimball models and 
the total electron number operator for the periodic An- 
derson model (we work in a canonical ensemble for the /- 
electrons in the Falicov-Kimball model, so no site-energy 
or chemical potential is needed for those particles). 

In Section II, we present a description of electronic 
charge reconstruction, which naturally occurs in any mul- 
tilayered device that can be grown for thermoelectric 
properties. This section briefly reviews the current sta- 
tus of such calculations, and describes their impact on 
the thermal transport; in particular, it fixes the notation 
for the internal electrostatic potentials associated with 
the electronic charge reconstruction. Section III provides 
the main arguments for developing the multilayered gen- 
eralization of the Jonson-Mahan theorem. Section IV 
applies the formalism to three classic experiments — the 
Peltier effect, the Seebeck effect, and the thermal conduc- 
tivity. Section V presents our conclusions and describes 
areas for further work. 
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II. ELECTRONIC CHARGE 
RECONSTRUCTION IN MULTILAYERED 
NANOSTRUCTURES 

The Schottky effect^, is a well-known effect in 
the semiconductor community, where charge is redis- 
tributed between a semiconductor and a metal at a 
semiconductor-metal interface due to a bulk chemical po- 
tential mismatch between the two materials. The charge 
rearrangement creates a screened dipole layer at the in- 
terface resulting in a final state with a static inhomoge- 
neous redistribution of charge through the system. Re- 
cently, the phenomenon has been revisited in the con- 
text of strongly correlated materials^i, where it has been 
called electronic charge reconstructioniSii^. If we imag- 
ine a multilayered nanostructure, composed of metallic 
leads sandwiching a barrier region, which is a strongly 
correlated material, then the chemical potential of the 
device is fixed by the chemical potential of the leads. 
If the chemical potential of the barrier is different, then 
the system must undergo an electronic charge reconstruc- 
tion. In particular, since the temperature dependence of 
the chemical potential should be different in the two dif- 
ferent materials, even if the chemical potentials match at 
one temperature, they will not match at other tempera- 
tures, and a charge redistribution will take place. 

In this work, we will focus on problems that have "mir- 
ror symmetry" for the leads, so the lead to the left is 
made of the same material as the lead to the right. In 
this case, the total Coulomb potential energy, due to all 
electric fields, goes to zero when one is far from the inter- 
faces, because all of the charge rearrangement is localized 
at the interfaces, and the whole system is charge neutral. 
If we were to examine systems with different materials for 
the left and right leads, then the electrochemical poten- 
tial of the system will be the average of the left and right 
bulk chemical potentials, which creates some additional 
complications, but does not change the basic strategies 
or formulas, although, the Seebeck effect needs to be de- 
fined and analyzed with carc^*. 

The approach to describe the electronic charge recon- 
struction is a semiclassical one. We solve the problem 
for local electron interactions exactly, but treat the long- 
range Coulomb interaction in a mean-field fashion. The 
strategy we use is to first calculate the electronic charge 
on each plane via an inhomogeneous Green's function ap- 
proach (in dynamical mean-field theory, the technique of 
Potthoff and Noltingi^ is used). If possible, one uses a 
Matsubara-frequency approach, because the numerics are 
usually under better control than real-axis approaches, 
but this is just a matter of convenience, not necessity. 
Next, we find the charge deviation on each plane; namely, 
we determine whether extra charge has entered or left the 
plane. Since the positive background charge of the ions 
remains the same, the charge deviation will give rise to 
an electric field. There are two different ways to treat 
this field. The simplest is to assume the electric charge 
is uniformly spread over the planeii. Then the electric 



field is constant, perpendicular to the plane, and pointing 
away from it in both directions if the net charge density 
is positive, while pointing toward the plane if the net 
charge density is negative. The second method uses the 
actual distribution of the ions, and the spatial profile of 
the electrons, if available, to calculate the charge^^. This 
approach is closer to an Ewald-like summation^^ of the 
charge densities. The two treatments should yield similar 
results. 

In this work, we will choose the "constant plane of 
charge" description for determining the electric fields. 
This allows us to determine simple analytic expressions 
for the electric fields — for example, the magnitude of the 
constant field, emanating from the ath plane of charge is 

where e < is the charge of the electron, pa is the 
quantum-mechanically calculated electron number den- 
sity at plane a, p^^^ is the bulk electron number density 
for the material that plane a is composed of (equal to 
the positive background charge on the plane), eq is the 
permittivity of free space, and e-^a is the relative per- 
mittivity of plane a. The contribution to the electric 
potential V'^{z) from this field satisfies 

E=-^V^{z). (7) 
dz 

Since the electric field is constant in magnitude, it 
is straightforward to compute the contribution to the 
Coulomb potential at plane (3 due to the change in the 
charge density at plane a (but one needs to keep track 
of the signs of the fields or equivalently the relative order 
of a with respect to (i): 

X <^ 0, /3 = a . (8) 

Note that if the relative permittivity is a constant, 
independent of the planes, then the potential energy is 
a linear function of the z-coordinate, proportional to 
— |Zq — z^l/ej. as one might expect. The reason why we 
need to sum over two terms in the summands in Eq. 
is because we envision the ath plane of charge to be in- 
finitesimally thick, and go through the lattice sites of 
plane a, but we assume the dielectric has a thickness of 
a and is centered around each plane of atoms. Hence, 
if the permittivity changes from one plane to another, 
a polarization charge develops halfway between the two 
planes where the dielectric is changing, and the electric 
field has a discontinuity at that point (i. e., at the posi- 
tion a + 1/2, see Fig. [TJ . 
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FIG. 1: Geometry taken for the classical electrostatics prob- 
lem. We show the blow up of two planes, a and a + 1. Assum- 
ing the net surface charge density on plane a is {pa— p^^^)a — 
Ua (located along the plane running through position a) and 
the relative permittivity is e^a (and similarly for the a -\- 1 
plane), then the change in polarization at the interface be- 
tween the two dielectric planes induces a polarization charge 
on the interface (denoted cTpoi) that leads to a discontinu- 
ous jump in the electric field halfway between the two lattice 
planes (at the position a + 1/2). Once the fields are known, 
we integrate to get the electric potentials. Note the disconti- 
nuity in the electric field occurs at the midpoint between the 
two lattice planes. 



It is actually the potential energy — |e|y^ = Va that 
shifts the chemical potential at each planar site. We de- 
fine a parameter 



eschot(a) 



2eoerQ 



(9) 



which controls how the extra charge density decays away 
from the interfaces. The parameter eschot has the miits 
of an energy multiplied by an area; the product of eschot 
with the local density of states has units of the inverse of 
a length, and this is what determines the decay length of 
the charge profile. Using this parameter, we can immedi- 
ately calculate the potential energy due to the Coulomb 
interaction (evaluated in a mean-field fashion) 



(10) 



T^j=a+i i[eschot(7) + eschot(7 - 1)]: /3 > a 
X <( 0, P = a 

J2j=a~i 5[eschot(7) + eschot(7 + 1)], P <a 

Note that a similar analysis can be carried out if one 
uses the Ewald-like technique for determining the charge 
reconstruction. 

These potential energies modify the Hamiltonian by 
the long-range Coulomb interaction of the charge recon- 
struction. The additional piece of the Hamiltonian (due 
to the charge rearrangement) is 



ct 2 (E plane 



(11) 



Hence, they can be treated by shifting the chemical po- 
tential — > /i — Vq, on each plane depending on what 
the Coulomb potential energy is for the given plane. For 
consistency, we must have that the potentials go to zero 
as we move far enough into either of the leads (for the 
mirror-symmetric case). This requirement enforces over- 
all charge conservation — any charge that moves out of 
the barrier remains in the leads, localized close to the in- 
terface, and vice versa. Of course, the potentials Va that 
appear in the electronic charge reconstruction Hamilto- 
nian in Eq. Ijlll) must be determined self-consistcntly. 
Achieving this goal requires care in setting up the itera- 
tive algorithm. 

There will be no electronic charge reconstruction if the 
chemical potentials in the bulk of both the leads and the 
barrier match. In order to have freedom to adjust the 
mismatch of the chemical potentials, we need to be able 
to change the value of the band zero of the barrier region 
relative to the band zero of the leads. This parameter is 
called AEpa, which vanishes in the leads, and is gener- 
ically a nonzero constant in the barrier (independent of 
the temperature or the charge rearrangement). Hence we 
add an additional term 



ffsct 



E E 

a 2 G plane 



(12) 



to the Hamiltonian. Although this term appears similar 
to the Coulomb potential term in Eq. Hll() , the key obser- 
vation is that this term is fixed and does not change with 
any parameters of the system, whereas the plane poten- 
tials Va need to be readjusted as the parameters change, 
to achieve a self-consistent solution of the problem. Note 
that we set AEpa = in the leads to the right and to 
the left. 

In this contribution, we will not discuss how to ac- 
tually solve for the electronic charge reconstruction in 
detail. One can imagine a number of different ap- 
proaches to this problem, ranging from direct means 
of solving the quantum-mechanical problem on finite- 
sized stacked planes, to other techniques like the inho- 
mogeneous dynamical mean-field theory approach. The 
DMFT approach has been quite successful in examining 
these kinds of problems, and the formalism only requires 
that the self-energy remain local (although it can vary 
from plane to plane). Then, by performing a Fourier 
transform to momentum space for the planar coordi- 
nates, one decouples the planar motion from the lon- 
gitudinal motion. Hence, the problem reduces to a se- 
ries of quasi-one-dimensional inhomogeneous problems, 
which can be solved by using the renormalized perturba- 
tion expansionii (sometimes called the quantum zipper 
algorithmic). Details for such an approach have already 
appearediiiC and are briefly reviewed below. 

The DMFT algorithm is given in Fig. |21 If there is a 
separate algorithm available for the Matsubara frequency 
Green's functions, then the upper left loop (which deter- 
mines the electronic charge reconstruction) is used on the 
imaginary axis, and we don't need it on the real axis. If 
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FIG. 2: Flow diagram for the DMFT algorithm in a multi- 
layered nanostructure with electronic charge reconstruction. 
We determine the charges on each plane, determine how they 
differ from the bulk charge on the plane to find the excess 
or deficit charge. Then we use classical electromagnetism to 
find the electric potentials on each plane and finally the con- 
tribution of the potential energy to the electrochemical po- 
tential on each plane. Then we average the potentials with 
a large damping factor so that the potentials are updated 
slowly. This is then input into the next loop of the main 
DMFT algorithm, which is unchanged from cases where there 
is no electronic charge reconstruction. 



culate the transport, we need to evaluate the real-axis 
results for the self-energies and Green's functions of a 
nanostructure with an electronic charge reconstruction. 
Unfortunately, the algorithms used when there is no elec- 
tronic charge reconstruction^'^' cannot be simply em- 
ployed for this case. The reason why is that the presence 
of the different potentials Va on each plane causes the 
nature of the integrands over the two-dimensional den- 
sity of states to have a different singular behavior than 
they had before. In a system without electronic charge 
reconstruction, the singularities in the integrand could be 
square-root-like, which are removed by a simple variable 
change using trigonometric or hyperbolic functions. Now, 
the singularities are poles (because the denominators are 
shifted by the potentials at a given plane, so they vanish 
at different energies, and give rise to a different singu- 
lar behavior), and we need to evaluate all integrals in a 
principal- value sense, where the real part is integrated 
with a symmetric grid around each pole, and the imag- 
inary part has a delta-function contribution that needs 
to be included. This is challenging to implement numer- 
ically, because the locations of the poles are different on 
different planes, and can vary from one iteration to the 
next. Details for how to deal with such a sophistication 
will appear elsewhere, since they are beyond the scope of 
this work. 



such an algorithm does not exist (such as when calcula- 
tions are performed with the numerical renormalization 
group), then we would use the entire loop on the real axis. 
The new steps to find the electronic charge reconstruc- 
tion are to first find the electron density on each plane. 
Then we subtract the bulk charge density of each plane 
to find the excess or deficit charge on the given plane. 
Once the change in charge density is known, we can cal- 
culate the electrical potential, and then the contribution 
to the potential energy. This gets added to the chemical 
potential to determine the electrochemical potential at 
each plane. 

When numerical results are generatedii*i^ (see those 
references for numerical issues in the algorithm), we find 
that usually the electronic charge reconstruction does not 
change significantly at low temperature, and that the size 
of the charge deviation grows as the mismatch between 
the chemical potentials grows (governed in part by the 
size of AEpa in the barrier). Such results are similar to 
what one would expect, but most of the calculations have 
taken place in systems where the charge density is not too 
sensitive to changes in the chemical potential. The effects 
may be different in systems with either Mott insulators, 
or doped Mott insulators, which can be brought close to 
the insulating phase via the electronic charge reconstruc- 
tion. In addition, electron-phonon coupled systems can 
develop a strong sensitivity of the charge to the chemical 
potential when the coupling is large, which may be an 
interesting case to examine as well. 

Ultimately, we are most interested in the transport of 
charge and heat through the device. In order to cal- 



III. PROOF OF THE INHOMOGENEOUS 
JONSON-MAHAN THEOREM 

It is important to examine how the linear-response 
transport formalismi^ is modified by the presence of 
an electronic charge reconstruction. We have taken the 
chemical potential as a constant throughout the multilay- 
ered nanostructure for thermodynamic equilibrium. One 
can directly show that the device carries no longitudi- 
nal charge current even though there are nonzero elec- 
tric fields arising from the electronic charge reconstruc- 
tion (see Appendix A for a proof when the self-energy 
is local). No current fiows because the putative current 
driven by the internal electric fields is canceled by an 
equal magnitude but oppositely directed current driven 
by the concentration gradients. The standard way to de- 
scribe this result is via a phenomcnological equation (for 
the case with no thermal gradicnts)^'^ 



if) = a^aapEp - a\e\^Vai3 

, 7 / "of 

e 



a X ^ 



(13) 



where Vap is the diffusion constant for Pick's law of dif- 
fusion^^-, and the second equality follows from the Ein- 
stein relation^? (or more correctly the Nernst-Einstein- 
Smoluchowski relation22iS4) which relates the diffusion 
constant to the conductivity via 

(Jap = e^Vapdpp/d^i] (14) 
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both quantities are matrices, with indices given by the 
planes of the multilayered nanostructure. The symbol 
= M ^ ^ct is called the electrochemical potential. The 
Einstein relation can be derived by relating the gradient 
with respect to the chemical potential to the gradient 
with respect to the number concentration via the chain 
rule: d^/dz = {dfi/dp)dp/dz, and the fact that the cur- 
rent vanishes in equilibrium.^'''. 

Eq. (|13ll implies that the condition for there to be no 
charge current is simply djl/dz — 0. The chemical po- 
tential is a constant, but it does vary with the filling, 
so if there is a change in electron concentration, then 
djl/dz = {dfi/dp)dp/dz ~ dV{z)/dz, so the force from 
the electric field will be balanced by the force from the 
change in electron concentration. In addition, note that 
the current vanishes no matter how large the variation 
in the concentration is (i. e., beyond the linear-response 
regime), so the conclusion is that the current generated 
by the internal electric field is always canceled by the 
current generated by the change in the electron concen- 
tration. Hence, for a linear-response treatment of trans- 
port, we can ignore the forces due to the internal electric 
fields and the concentration gradients, because they al- 
ways cancel, and we can limit our focus to the effects of 
the external electric field only. This then implies that 
all of the analysis performed previously for the charge 
currenbi^ continues to be valid, and because the form of 
the charge current is unchanged when we have electronic 
charge reconstruction, the Kubo formula is identical as it 
was before (with the effects of the potentials Va included, 
of course). 

The basic observation needed for a thermoelectric de- 
vice is that there is a difference between the weighting 
factors that determine the bulk charge current and heat 
current. The charge current is weighted by the electron 
velocity, while the heat current is weighted by the velocity 
multiplied by the kinetic energy minus the chemical po- 
tential plus a term from the potential energy. Hence, one 
can create charge current without heat current, or vice 
versa; by carefully engineering the way electrons move 
through the device, one can control both the energy and 
charge flow, which is useful for different types of appli- 
cations like refrigeration or power generation. A typical 
device has two legs, one using electrons as the charge car- 
riers and one using holes as the charge carriers. Current 
flows through the device in a loop, but the net heat flows 
in one direction only, which allows the device to function. 

In this contribution, we concentrate on multilayered 
nanostructures, which can be used to compose one of 
the legs of the thermoelectric device. We also concen- 
trate solely on electronic transport mechanisms. In most 
thermoelectrics, the thermal conductivity from phonons 
can be large enough to significantly reduce the figure-of- 
merit. It is expected that the phonon thermal conduc- 
tivity will be further reduced in a nanostructure, because 
the interfaces in the nanostructures will cause significant 
phonon scattering if the masses of the ions in the different 
materials have a large mismatch?;, but we do not discuss 



this issue further. 

There is no simple way to derive the response of a 
strongly correlated system to both electrical fields and 
thermal gradients. The reason why is that the thermal 
gradient cannot be added as a field to the Hamiltonian 
like the electric field can, hence there is no way to fol- 
low the simple Kubo response theory developed for the 
charge current in an electric field (because the linear- 
response approach evaluates correlation functions at a 
fixed temperature, and a variation of the temperature 
with position is problematic to include within the for- 
malism). Luttinger sorted out a reasonable plan of ac- 
tion for how one can nevertheless proceedS^. We couple 
a fictitious field to the heat-current operator, analogous 
to the vector potential that couples to the charge current 
operator, and determine the linear response with respect 
to both fields. Then, we compare the Kubo response to a 
phenomenological set of equations that relate the charge 
and heat currents to the electric field and the gradient 
of the temperature. We then identify the relevant trans- 
port coefficients and how they are expressed in terms of 
correlation functions. 

In multilayered nanostructures, there always is an elec- 
tronic charge reconstruction, because the bulk chemical 
potentials for the leads and the barrier will have differ- 
ent T dependence, and hence cannot always be equal 
(the only exception is for particle-hole symmetry at half- 
filling, but there the thermopower vanishes, so that case 
is uninteresting for thermal transport). Hence the Hamil- 
tonian must be modified to include the potential energy 
Va on each plane, and the band offsets AEp, as described 
in Sec. II [i. e., we add Y^aii^a - AEFa)cliCa, to H]. 
The band offsets are independent of T, and represent 
the difference in the band zeroes for the leads and the 
material placed at plane a. The potential energies Va 
do depend on T, but they do not create any currents, 
because they correspond to the static potential associ- 
ated with the electronic charge reconstruction (and the 
diffusion current generated by the change in electron con- 
centration cancels the current from the internal electric 
field; see Appendix A). But the Coulomb potentials do 
create internal electric fields that maintain the electronic 
charge redistribution amongst the planes. 

The phenomenological study of currents caused by ex- 
ternal electric fields or temperature gradients has been 
examined since the early 1800s. It was found that an 
electric field can drive a charge current (which is essen- 
tially Ohm's law2& with the conductivity as the phe- 
nomenological constant) and it can drive a heat cur- 
rent because the electrons carry heat with them as they 
move through the material (this phenomenon is called the 
Peltier effect^^). Similarly, a temperature gradient can 
drive heat conduction with the phenomenological ther- 
mal conductivity (called Fourier's law^^), and because 
the electronic contribution to the heat current generi- 
cally carries charge, a temperature gradient can gener- 
ate a charge current (called the Seebeck effecliSS). The 
phenomenological equations for the (linear response) Ion- 
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gitudinal transport in a multilayered nanostructure are 
then (Jq, is the longitudinal number current, is the lon- 
gitudinal charge current, and is the longitudinal heat 
current) 



e a 



e|(ja) 



dT a 



|e|a^ Li2q/3 



~ Tf} 

dT a 



(15) 



Ie|£;fi 



(16) 

where the indices a and /3 denote the planar sites (or 
the midpoint between planar sites, as clarified below), 
the term (T/3+1 — Tfj)/a is the discretized approximation 
to the temperature gradient and the L^j coefficients can 
be thought of as the phenomenological parameters. We 
define the symbol /i^g = /i — Vg + AEpfj, which may be 
thought of as the "local chemical potential" for plane 
j3. The origin of the temperature derivative of enter- 
ing into the phenomenological equations arises from the 
conventional V/i term, which becomes VTdfi/dT when 
the system is placed in a thermal gradient. The spatial 
derivative of the terms does not drive any current, 
because it cancels with the current driven by the equi- 
librium concentration gradient (which we did not include 
in the above phenomenological equations) , so the electric 
field Ef3 is the external field applied to the device (this 
is valid only in the linear-response regime of a small ex- 
ternal electric field). Note, that there is a simple way to 
understand the signs that appear in Eqs. H15|) and (|16|l . 
First consider the external electric field, which can be 
written as the negative gradient of the electric potential. 
The current (whether of electrons or of holes), always 
runs down the potential hill. Since the conductivity is 
always positive, the first term in Eq. H15|l must have a 
positive sign. The thermoelectric number current also 
runs downhill, so it is proportional to the negative tem- 
perature gradient. For electrons, the charge current is 
— |e| times the number current, which gives rise to the 
positive sign for the last term in Eq. (|15() . Similarly, the 
thermal conductivity runs down the temperature "hill" , 
so it has a negative sign in front of it. The Peltier effect 
term is the hardest to understand, but because the elec- 
trons are negatively charged, they actually move up the 
potential hill (the charge current runs down the hill be- 
cause the electrons are negatively charged), so the heat 
is carried up the hill, and hence there is a minus sign in 
front of the term (recall the electric field is the negative 
gradient of the potential). 

Our next step is to determine how to represent the 
thermal transport coefficients Lij in terms of many-body 



correlation functions. This has already been done for 
the first coefficient-'^^, which is proportional to the con- 
ductivity matrix, and is represented by a current-current 
correlation function: aa/3 = c^Lhq^ (the modification of 
the Hamiltonian by the electronic charge reconstruction 
has no effect on the form of the charge current, or on the 
form of the correlations functions, but obviously creates 
additional scattering). Since this coefficient arises from 
an electric field, which can be added to the Hamilto- 
nian, the derivation is rigorous. Similarly, if we follow all 
the steps in Ref. |3 that led up to the derivation of the 
conductivity matrix, but we examined the expectation 
value of the heat-current operator instead of the charge- 
current operator, we would find that the L21 correlation 
function was identical to the Ln correlation function ex- 
cept that it is a heat-current-charge-current correlation 
function instead of a charge-current-charge-current cor- 
relation function. 

As we discussed above, there is no complete the- 
ory to determine the L12 and L22 coefficients for the 
phenomenological transport equations. But, classical 
noncquilibrium statistical mechanics has proved that 
there is a reciprocal relation between the "cross" terms 
in the transport equations^S.. Written in the form we 
have them, this relation says that L21 = Li2. Know- 
ing the form for L21, we then conclude that L12 is the 
charge-current-heat-current correlation function. Keep- 
ing within this same vein, the natural conclusion is that 
the final transport coefficient L22 is a heat-current-heat- 
current correlation function (but there is no rigorous 
derivation of this result) . 

In order to derive the local charge and heat current 
operators, we must formulate the transport problem in 
real space. Unlike the bulk case, where the procedure is 
completely well-defined, there are a number of different 
possible ways to try to derive the local current opera- 
tors. The bulk number current operator is found by tak- 
ing the commutator of the number polarization operator 
with the Hamiltonian; this guarantees that the equation 
of continuity holds, and it also implies that the number 
current is conserved through the system. The number 
polarization operator is 



TTimmbor — \^ JJ .(J r + f'' f ■) 



(17) 



where we dropped the spin index for simplicity, and 
where Rai is the position vector of the site labeled by 
ai. The p f term enters for the Falicov-Kimball and pe- 
riodic Anderson models, but not for the Hubbard model. 
Since the number polarization operator depends only on 
the number operators, it commutes with all number op- 
erators in the interaction Hamiltonian. It turns out that 
the form in Eq. I|17|) also commutes with the hybridiza- 
tion term in the periodic Anderson model because the 
hybridization is on-site only. Hence, the bulk number 
current operator is the same for all models we examine 
here, and arises solely from the commutator of the z- 
component of the polarization operator with the hopping 
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Hamiltonian. Performing the commutator is straightfor- 
ward, and leads to j = j^, with 

iG plane 

Note that the subscript a on the current operator de- 
notes the total current operator flowing through the ath 
plane, and does not indicate a Cartesian coordinate of the 
current operator; the current operator is always taken in 
the z-direction for the longitudinal flow. The current at 
plane a is thus defined to be the total number of elec- 
trons flowing to the left minus the total flowing to the 
right (here the current operator at plane a is determined 
by the number of electrons flowing to the right or to the 
left through the ath and a + 1st planes). 

A comment is in order about the choice given in 
Eq. 1)18(1 for the current associated with the ath plane. 
Note that the form chosen is not the same as the choice 
that would arise from taking the commutator of the "lo- 
cal" polarization operator (at the ath plane) with the 
Hamiltonian. The direct result from the commutator 

ja = X/iGplanc ^o^{c^^C^^ + faifai)] 



ja — i [^aa+l('^I+l 
iGplanc 



ai'-'a+li) 

+ ^a-la{^a-li'^ai ~ '^ai'^a-li)\^on (19) 



does not seem reasonable, because it is weighted by the 
z-coordinate of the ath plane, rather than involving the 
difference of currents moving in opposite directions (at 
the ath plane). When we have full translational symme- 
try, we derive the conventional form for the current oper- 
ator by shifting the spatial index of one of the terms, to 
explicitly carry out the cancellation of the spatial coordi- 
nates (just take the summation of the above result over 
a, and shift a a+1 in the last two terms). More reflec- 
tion on this issue, shows that the explicit form of the local 
current operator that enters the Kubo formula actually 
originates from the coupling term — j • A term that corre- 
sponds to the perturbation of the Hamiltonian due to the 
electric field in a gauge where the scalar potential van- 
ishes; this is because we evaluate the expectation value 
of the total current with the perturbation of the Hamil- 
tonian due to the external field and that field enters via 
the vector potential value at a specific plane. Hence the 
conductivity matrix is defined from the piece of the total 
current operator that couples to the field at plane a and, 
since the total current will be the sum of the currents at 
each plane, the current-current correlation function for 
the conductivity matrix involves the local current opera- 
tors that couple to the vector potential. Thus, we choose 
the perturbation of the Hamiltonian to be 

V-'it) = -i|e|a^t^„+i(cJ,+iiC„, -c^c„+n)^a(t), 

(20) 

where we have taken the vector potential along the z- 
direction, and independent of the intraplane coordinates. 



because the field is uniform for each plane. We feel this 
choice makes good physical sense because we couple the 
vector potential to the physical current between the ath 
and a 4- 1st planes. Alternatively, one can view this as a 
coupling of the current between the ath and a -I- 1st plane 
to the electric vector potential located halfway between 
those two planes (in this interpretation, we would use 
[Aa + Aa+i\/2 as the coupling field). Finally, one can 
take a symmetrized version of the local current operator 
to be 

iG plane 

— iatj^^^i ^ {^'^a + li^ai ~ cL'^Q+lj)/^: (21) 
iG plane 

corresponding to the average of the currents located just 
to the left and to the right of plane a. This choice sounds 
like the most physical choice, but the calculations for 
it are somewhat more complicated, and it is not likely 
the end results are too different from our first choice. 
The difference between the two choices is actually quite 
simple. In the first approach, one should envision the 
spatial indices a and /3 to correspond to Za + a/2 and 
z/3 -I- a/2; that is, they are shifted to the right by half the 
distance between the planes. In the second, symmetrized 
approach, the a and [3 indices denote the planar indices. 
For this reason, we don't expect the final results to be too 
different for either approach. For simplicity, we choose to 
take the current operator to be the current between the 
ath and a -I- 1st planes for our derivations below, and we 
discuss how to get the corresponding symmetrized results 
at the end. 

The calculation of the local heat current operator is 
more complicated. We adopt the same strategy as be- 
fore though — first calculate the bulk operator, and then 
extract a reasonable choice for the local operator. To 
calculate the bulk heat current operator, we first need 
to determine the energy polarization operator. This is 
similar to the number polarization operator, except it is 
weighted by the piece of the Hamiltonian associated with 
the ai position. This is easy to do for the interaction 
and hybridization terms, the Coulomb potential energy 
terms, and the band offset terms, which are local, but 
is complicated for the hopping terms, which involve two 
lattice sites. The procedure that is used is to associate 
half of the hopping term between the two sites with the 
local Hamiltonian at each of those two sites. The energy 
polarization term then becomes 



H^ 



E E 

a ■iGplanc 
^charge ai 



2 5^ 5^ ^hop aif3j 
(3 plane 



n 



offset ai_ 



Rn 



(22) 



where the hopping piece is divided into two as described 
above, and the interaction piece includes all the local 
parts of the interacting Hamiltonian associated with each 
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lattice site. The bulk energy current operator is = groups the terms together to find how to make an edu- 

i[H, n^], and the heat-current operator is = — fij cated guess for the local heat current operator. The final 

because the heat is the energy measured relative to the results that wo have are summarized below. These are 

chemical potential. The commutator is tedious to work the proper local heat current operators needed to satisfy 

out, but just involves straightforward algebra. When it is the Jonson-Mahan theorem, as described below. In all 

finished, we have an expression for the bulk heat current, cases, we have j'^ = j^. Note that we can form the 

which can be organized into summations that involve a symmetric version of the heat current operator as well, if 

plane a and the plane to the right (there is also a hopping desired, but it is even more complex, 
term involving operators at the a + 2 plane). One simply 



For the Hubbard model, we have 



E 

iG plane, (7 



iG plane, (7 

1 

2 



2''ce— la / V \^OL-\-li<7^oc—li<7 *^a — licr^^a+lifr) 
iG plane, (7 



iGplane,(T ) 

where a = —a denotes the spin state opposite to a. For the Falicov-Kimball model, we find 

V ijGplane iGplane 



(23) 



E 

iG plane 



(24) 



For the periodic Anderson model the commutation of the Hamiltonian with the energy polarization operator 

2^a+la+2 E^ (^a+2i<7Cajo- — C^ai<T'^a+2ii7) ~ 2^a-la E^ i^a+lia'^a-lia " 'it-lia'^a+licr) 



gives 



3a = iat, 



+ E 

iG plane, cr 



iG plane, cr 

1 

2 



iG plane, cr 



+ iaV^^ 2 y ^ (/ai(T''a + lio- ''a+licr/aio-) "I" (/ai(T''a-licr '^a-licrfaia) 

iG plane, cr 



iG plane, cr 



(25) 



The heat current operator depends on the model being tract the chemical potential multiplied by the number 

examined, because it involves commutators of the poten- current from the energy current to get the heat current, 
tial energy with the energy polarization. We also sub- One might have thought we should subtract the "local 
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chemical potential" multiplied by the local number cur- 
rent operator, but that would remove the extra terms in 
the heat current arising from the electronic charge re- 
construction; one could have grouped those terms into 
either the Hamiltonian or the local chemical potential — 
we chose the former, so we subtract only fij. 

Now we need to determine the dc limit of the cor- 
relation functions on the real axis. The analytic- 
continuation procedure is identical to that for the bulk 
case. We start by defining a polarization operator on the 
imaginary axis, then we analytically continue to the real 
axis, we form the relevant transport coefficient, and then 
we take the limit of the frequency going to zero. We de- 
note four polarization operators by Lijapiivi) according 
to 



and the transport coefficients satisfy 



Lijap = Imi — [Ly a/3 (i^ + «0+) - 0/3(1/ + iO )] 



lim Re[-iLyQ^(z/)/z/] 



(27) 



Lllal3{m) = 



dre-'-(r.j„(r)j^(0)), 
'dre-'-(T.iJr)jQ(0)), 
'dre-'-(r.jQ(r)j^(0)), 
'dre-'^(r.jQ(r)jQ(0)}, (26) 



(the ij subscripts here are 1 or 2, and not the planar 
site indices). The generic notation 0{t) = exp[(?i — 
IJ,N)t]0 e.yi\)[—{H — fiM)T] is used to indicate the time 
dependence of the operators in Eq. I|26|l . The Jonson- 
Mahan theorem^''^ can be straightforwardly generalized 
to treat this case. Begin by defining a generalized func- 
tion 



-Fa/3(Ti,r2,r3,r4) = ^T^iat^„+i ^ [cl+ii(ri)c„i(r2) - c^(Ti)c„+i,(r2 

X iatj^P+l [4+lj(^3)%(^4) -4j(^3)C;3+1j(t4) )■ 

j'G plane 

Next, we determine the polarization operators by taking the appropriate hmits and derivatives. Namely, 



(28) 
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Ll2ap{il^l) = 


/ drie 




Jo 






L21al3{m) = 


1 dTie 




Jo 






L22af}{m) = 


/ dTie 




Jo 
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9t2 




( d 


d 






dT2 



T3=Q,Ti=0- 



^'q/3(ti,T2,0,0 



^q/3(ti, T2, T3, r4) 



(29) 



This result holds because the [dr — 9t')/2 operator converts the local number current operator into the local heat 
current operator. To see this, we simply compute 



iGplanc 

= *<a+l {[^-Ai-A/',cl + nW]c«.(T)+cl+i,(T)[H-/iA/-,C„(T)]}, 

z£ plane 



(30) 



which can be shown to be equal to when the commuta- tors are evaluated. It is this critical identity that connects 
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the local number and heat current operators that is a re- 
quirement for the formalism to satisfy the Jonson-Mahan 
theorem. The analytic continuation is complex, because 
it involves four-time functions in the general case, and 
a detailed proof of the Jonson-Mahan theorem appears 
in Appendix B. Instead, we provide a direct construc- 
tive proof in DMFT here, where we neglect the vertex 
corrections. This is just a heuristic approach to the full 
problem. 

The first step is to evaluate the expectation values of 
the Fermionic operators (in the definition of F) via con- 
tractions, because we neglect the vertex corrections. This 
yields 

{Gfja+lji{TA ~ Ti)Gal3+li3{T2 — T3) 

ijG plane 

— Gp+la + lji{T4 — Ti)Gal3ij{T2 — T3) 

— Gpaji{T4 — Ti)Ga+ip+lij{T2 — T^) 

+ Gi3+iaji{T4 — Ti)Ga+lf3ij{T2 - T3)} . (31) 

Next, we need to determine a spectral representation for 
the off-diagonal Green's function. Using the fact that 



Gal3ij{z) = J duj 



lmGaf3ij{i^) 



(32) 



with z in the upper half plane (which can be shown by 
using the Lehmann representation), says that 



GafUjir) = J dcjTy^^ 



IUJ„ 



-lu\Gapij{uj). (33) 



Now we convert the sum over Matsubara frequencies into 
a contour integral (that surrounds each Matsubara fre- 
quency, but does not cross the real axis — the contour is 
then deformed into two contours, one running just above 
and the other just below the real axis), but we must be 
careful to ensure that the procedure is well-defined. If 
T < 0, then 



T 



i 


/ dz — 

Jc 


ZT 


i 








/ dze 
'J —00 


-^V(^) 




1 


1 


z + 




z — iO+ — Lu 



7(c^). 



(34) 



This result is well-defined because the Fermi factor pro- 
vides convergence (asymptotically like exp[— /?z]) for z — > 
00 and the exp[— zr] term provides boundedness for 
z — > —00 when r < 0. Since 1 — f{z) has the same 
poles as f{z) on the imaginary axis, with residues that 
have the opposite sign, and it behaves like exp[/3z] for 
z — oo, one finds 



e---[l-/(c.)]. 



(35) 



for r > 0. The results in Eqs. (|34() and (|35|l can then be 
substituted into Eq. H33|l to get the final formula for the 
off-diagonal Green's function 



GaPij (t ) 



-i / dujlxnGo,p,j{Lu)e-'^-[l - f{io)], T > 
--/ dujl-mGai3i]{^^y 



)e -^[-/(c^)], T<0. 

(36) 

Now we note that we can restrict ourselves to the case 
Ti > T2 > Ta > T4 without loss of generality, because 
that is the ordering needed to get the relevant correlation 
functions. Then we employ Eq. (|36|l in Eq. (|31|l and 
use the fact that the summations over the spatial indices 
for the planes can be Fourier transformed, and then the 
momentum summation can be replaced by an integration 
over the two-dimensional density of states, to yield 



a j_ 



(37) 



duj J duj' J rfellp^''(e") 
X f{uj)[l - /(w')]e-'^(^''-^i)-'^'(^=-^3) 
X |lmG/3a(ell , Lu)lmGa+if}+i (e" , uj') 

+ImG/3+ia+i(e" , t^)ImGQ/3(e" , w') 
~lraG/3a+i (e" , t^)ImGa/3+i (e" , uj') 
-ImG/3+iQ (e" , w)ImGQ+i/3(e" , uj') | 

Now we can evaluate the polarizations, and directly per- 
form the analytic continuation. We Fourier transform the 
expression in Eq. (|37|l to get the Matsubara frequency 
representation. Then we replace ii^i hy v + iO^, we con- 
struct the transport coefficients on the real axis, and we 
finally take the limit ^ to get the dc response. The 
factor {dr — dr')/2 gives a factor of {uj + uj')/2 which goes 
to (w-f J//2) after integrating over the delta function that 
arises in the analytic continuation. Setting v — gives 
an extra power of uj in the integrand for each derivative 
factor in the response coefficient. The end result is 



± 



± 



duj — 



dm 

dio 



(38) 



dellp2'^(ell) 

{ lmGpaieKuj)lmGa+i(3+i{e\uj) 

+ ImGQ/5(e",Lj)ImG/3+iQ+i(e",a;) 

- ImG^jQ+i(ell, w)ImGQ/5+i(e",t^) 

- ImGQ+i/3(e",Ci;)ImG/3+ia(e",t^)}- 

This is the generalized Jonson-Mahan theorem for in- 
homogeneous systems described by DMFT with vertex 
corrections neglected. Note that the equality of L12 with 
L21 is the Onsager reciprocal relationSS,. 

If one wants to work with symmetrized currents rather 
than the currents between the a and a + 1st planes, then 
the Kubo formulas will be changed slightly to take into 
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account the symmetrized current operators. These can 
be constructed directly from the correlation functions al- 
ready illustrated above, and it is a simple exercise to take 
care of the relevant bookkeeping; we leave such details to 
the reader. 



IV. ANALYSIS OF EXPERIMENTS 

With the expressions for the phenomenological coef- 
ficients that appear in Eqs. IjlSI) and Hlt)|) determined, 
we now can move onto evaluating the transport in dif- 
ferent cases of interest. The first point that needs to 
be emphasized is that the total number of electrons is 
always conserved in the system, so the charge current 
is conserved, and cannot change from plane to plane 
(Jq) = {jp)- There is no such conservation law for the 
heat current though, because the electrons can change 
the amount of heat that they carry depending on their 
local environment. Hence, it is the boundary conditions 
that we impose upon the heat current that determines 
how it behaves in a multilayered nanostructure. This 
point will become important as we analyze different ex- 
perimental situations. 




iQL-^«-^ Iqr-^ 

I I I I ♦ I 

AIq, ^\o2 ^Jq3 A'Q4 AjQ5 

FIG. 3: Schematic diagram of the heat transfers in the Peltier 
effect. The dots refer to different planes. A heat current is 
incident from the left. As we go from one plane to another, 
the heat current changes, as heat is transfered to or from a 
reservoir to maintain the system at a constant temperature. 
For example, we can examine the total heat current transfered 
to the reservoirs {jqr — Jql), or we can examine the average 
heat current that flows through the device "^Za JQ'^/-^ ' where 
A'^ is the number of planes involved in the heat transfer. 

The first experiment we would like to analyze is the 
Peltier effect in a multilayered nanostructure. We imag- 
ine that the nanostructure is attached to a bath that 
maintains the entire structure at a fixed temperature, 
and we then turn on an external electric field. The Peltier 
effect is the ratio of the heat current to the charge cur- 
rent. A moment's reflection will show that the heat cur- 
rent is not necessarily conserved in this system, because 
we have to exchange heat with the reservoir to maintain 
a constant temperature profile. Hence, it isn't even ob- 
vious what ratio should be taken for the Peltier effect — 
the average heat current for the entire device over the 
charge current, the total change in the heat current over 
the charge current, or the heat current transfered to the 



reservoir over the charge current. We now show how to 
determine all three of these results. 

The starting point is the transport equations [ (|15|l and 
(|16|l ] with Ta = T independent of the plane number. We 
first determine the electric field by multiplying both sides 
of Eq. (|15|l by the inverse Ln matrix. Since the charge 
current is independent of the plane index, we find the 
electric field satisfies 

^" = iEM)a,(/)- (39) 

Integrating the electric field over the z-coordinate, then 
yields the voltage across the device, which allows us to 
extract the resistance-unit-cell-area product via Ohm's 
law 

«"«' = ^E(^n)a/5- (40) 

To find the heat current, we substitute the value of the 
electric field into Eq. Ijltill . which yields 

(j-) = -oE^2i./5(ir/);3,(/)- (41) 

This is all we need to analyze the Peltier effect of a nanos- 
tructure. Note that the heat current generically will have 
a dependence, and hence will vary from plane to plane 
(see Fig. El. 

The first question we can ask is how much heat is lost 
or gained by the reservoir that is attached to the device 
to maintain isothermal conditions. This is determined by 
the ratio of the difference in the heat current at the right 
and the heat current at the left to the charge current. In 
equations, 

A(jQ> ^ (.7g) - (.7g) 
if) if) 

- -oE(^2ifl/3-W/3)(iri')^,- (42) 

This would measure the net cooling or heating of the 
reservoir by the device as the charge current flows. Sim- 
ilarly, we could measure the average heat flow carried 
through the device 

%f = ^O^E^21a/3(iri')^^, (43) 

where N is the number of terms taken in the summation 
over the index a. This expression is analogous to the bulk 
Peltier effect, which measures the ratio of the heat to 
charge current flows (which are independent of position 
in a bulk system in linear response). 

Next we examine the Seebeck effect and a thermal con- 
ductivity experiment. In both cases we work with an 
open circuit, so the total charge current vanishes (j'^) = 0. 
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drop across the device. We also need the temperature 
profile. Substituting Eq. (g^ into Eq. (fTC|) . and noting 
that the heat current is conserved if the device is isolated 
and in the steady state (implying heat cannot be trans- 
fcrcd out of any plane, recall that the Joule heating is a 
nonlinear effect) because the system develops a tempera- 
ture profile so that the heat current is conserved through 
the device. In this case, we can evaluate the temperature 
profile, which satisfies 



a+l 



-tJ2{m- 



lap 



(45) 



FIG. 4: Schematic diagram for how to measure the (relative) 
Seebeck effect. The metallic leads are composed of the same 
material, and a voltage probe is placed across the two ends 
which both are fixed at a temperature of Tq. The voltage 
across those probes is Vo- Since heat current will flow across 
the voltmeter if both ends are not at the same temperature, 
there is no way to directly measure the desired voltage V. 
But, since the change in voltage in the metallic lead in going 
from To to T on the left hand side is exactly canceled by the 
change in voltage in going from T to To on the right hand side, 
we find the difference between the voltage V and Vo is just 
equal to the Seebeck coefflcient of the metallic lead multiplied 
by AT. Hence, since a measurement uses Vo instead of V, the 
Seebeck coefficient of the barrier is measured relative to the 
Seebeck coefficient of the metallic lead. 



The Seebeck measurement is subtle, because we don't 
want to measure the voltage difference with probes at dif- 
ferent temperatures, because there will be a contribution 
from the VTdji/dT terms to the voltage drop (and there 
may be a thermal link allowing heat to flow through the 
voltage probe). An actual experiment uses thermocouple 
probes, where one end of the probe is placed on the sam- 
ple, and the other is placed in a constant Tq bath. Two 
probes are needed to measure the voltage change and the 
temperature at two points along the sample. The net 
thermopower is measured relative to the thermopower of 
the metal used in one of the legs of the thermo coup le 
(typically copper). For details, see Refs. HO and l3lt a 
simpler schematic picture of this issue is shown in Fig. 21 
Alternatively, we can imagine the lead to the left placed 
in a bath at temperature Tg, the interface plane on the 
left held at temperature T, the interface on the right held 
at temperature T + AT, and the lead to the right held 
at temperature Tq. The net effect on our analysis, if we 
assume the thermopower of copper can be neglected (or 
of the ballistic lead in the alternative picture) , is that we 
neglect the dfj,a/dT terms in our analysis (because the 
chemical potential at the probes is at a constant temper- 
ature when the potential difference is measured). With 
these caveats in mind, using Eq. (|15|l . we find 



^ E ii2^,(r^+i - T,). (44) 



alelT 
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with the matrix M defined to be 

Ma(l = L22al3 ^ ^ -^21a7 (L^i) g L 



(46) 



7(5 



Now we can sum Eq. (|45ll over a to get the tempera- 
ture difference over the device. Hence the Seebeck effect 
becomes 



Ay 1 Y.aP'^5iLii)al3Li2p^M^l 

\e\T Eaf^Mj 



(47) 



Note that this is not equal to 1/T times the Peltier coef- 
ficient as in the bulk [see Eqs. (|42|l and (|43|) ]. Instead, we 
have a weighting of the L12 to Ln ratio by the matrix Af , 
which is related to the thermal conductivity. This factor 
cancels in the bulk, where the M matrix depends on the 
difference of the spatial coordinates, and the q = re- 
sponse is independent of M because the common factor 
in the Fourier transform will cancel out (as can easily be 
proved by invoking the convolution theorem). If we do 
not measure AV via thermocouples at constant T, then 
the AV term is modified by a contribution from dfia/dT. 
Wc do not discuss that modification here, because it is 
not normally a technique used in measurements^^. 

The thermal conductance is evaluated in a similar way, 
but does not require any subtlety in the measurement. 
We also work in an open circuit, and the heat current 
is conserved, because we isolate the system. Now we 
measure the ratio of the heat current to the temperature 
difference to find that the thermal conductance per unit 
area K satisfies 



K = - 



1 



AT tj:^^{m-^w 

and the thermal resistance-area product becomes 
i?,ha^ = T^(M-i)^^. 

a/3 



(48) 



(49) 



Multiplying by a and summing over a yields the voltage 



Given all of the phenomenological parameters that en- 
ter into the transport of a nanostructure, we are now in 
a position to be able to evaluate things like the efficiency 
of a refrigerator, or of a power generator (this requires 
evaluating heat flow while charge current is flowing, and 
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is more complex than the cases considered here) . The fi- 
nal equations that result are quite complicated, and will 
not be shown here. Note that it is necessary to perform 
such an exercise here, because in these nanostructure de- 
vices, the thermoelectric cooling or power generation is 
not determined solely by the bulk figure of merit of the 
constituent pieces. When quantum effects enter due to 
nanoscale structures, the situation is more complicated. 
But one can define an effective figure-of-merit by con- 
structing an effective Lorenz number from the ratio of 
the thermal to the charge resistance and then evaluate 
an effective figure-of-merit from the Seebeck coefficient 
and the effective Lorenz number. Even if this is done, it 
is not the same as the calculation of the efficiency of a 
real device, which simply is more complicated. 



CONCLUSIONS 



In this work we have shown how to derive the for- 
malism for evaluating the electronic contribution to the 
charge and thermal transport of a strongly correlated 
nanostructure by determining the appropriate Kubo for- 
mulas for the transport coefficients. This analysis re- 
quires us to properly determine the local current opera- 
tors, which we do via a heuristic argument. Using these 
specific local current operators allows us to show that the 
Kubo formulas for the various heat transport coefficients 
are related by a generalization of the Jonson-Mahan the- 
orem to inhomogeneous systems. We also describe how 
nanostructures that will be used for heat transport must 
have an associated electronic charge reconstruction, and 
we sketched how to solve for that reconstruction using 
the DMFT approximation in the nanostructure. We il- 
lustrated our results for models described by the Hub- 
bard model, the Falicov-Kimball model, and the periodic 
Anderson model, but they should hold for any model that 
involves only local interactions. 

It will be interesting to now solve some numerical prob- 
lems and investigate the thermal transport coefficients for 
these systems. Devices of particular interest are those 
that have strongly correlated electrons, such as systems 
with a Mott insulator, a doped Mott insulator, or Kondo 
metals in them. Solving such problems is easiest to do 
for the Falicov-Kimball model because it is the easiest 
model to solve numerically, but it should be feasible to 
investigate the Hubbard model and the periodic Ander- 
son model as well using numerical renormalization group 
techniques, or quantum Monte Carlo plus maximum en- 
tropy analytic continuations. Such work will be presented 
elsewhere. 



APPENDIX A: VANISHING OF CURRENTS FOR 
AN EQUILIBRIUM ELECTRONIC CHARGE 
RECONSTRUCTION 

We prove that the charge and heat current expectation 
values vanish in equilibrium when there is an electronic 
charge reconstruction, in the case where the self-energy is 
local. The result should also hold for the nonlocal case, 
by purely physical reasons, but we are not aware of a 
simple way to prove this result in the general case. 

Because the multilayered system is translationally in- 
variant in each plane, we can use a mixed basis, where we 
have a two-dimensional momentum describing the planar 
degrees of freedom, and we use real-space to describe the 
inhomogeneous z-direction^^. Then, if we assume the 
self-energy is a local function Sq,, that can vary from 
one plane to another, we have the local Green's function 
satisfiesi^ 

G„„(kll,o;) = l/{i„(kll,c^)-f i?„(k'l,o.) (Al) 

where is the two-dimensional (planar) bandstructure 
on plane a, and the left La and right Ra functions are 
defined below. The local Green's function on each plane 
is then found by summing over the two-dimensional mo- 
menta, which can be replaced by an integral over the 
two-dimensional density of states: 

GaaH = j deip^\ei)G^a{ei,u:), (A2) 

The left function is defined to be 

La_„(kii,w) = — -t (A3) 

GaQ_„(kll,tj) 

and it satisfies the recurrence relation 

ia_„(kll,cj) = + /i - + A£;Fa - Sa_„(w) - e|[_^j^| 

a — '/la — 71— 1 a—n—la — n 



LQ_„_i(kll,a;) 



(A4) 



We solve the recurrence relation by starting with the re- 
sult for i_oo, and then iterating Eq. (|A4I) . In a similar 
fashion, we define a right function and a recurrence rela- 
tion to the right, with the right function satisfying 



i?a+„(kll, W 



Gc(Q+n-l(k" , ijj)ta+n-la 



(A5) 



Gaa+n(kll, w) 

and the recurrence relation being 

i?a+n(kll,Cj) = W + /i - Va + AE'pa - SQ+,i(t^) - e|[^ 



t 



a+na+n-t-l a-t-n-t-la+n 



-Ra+ji+l(kll , 



(A6) 



We solve the right recurrence relation by starting with 
the result for i?oo, and then iterating Eq. (jA6p . 
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In order to determine the current, we need to exam- 
ine the ofF-diagonal, nearest neighbor Green's functions, 
which satisfy 

Ga+iQ^k", , , (A7) 

La(kll,cj) 

and 

Gc«a+i (kii , iS) = — — . (A8) 

i?Q+i(kll,a;) 

Using the recursion relations for the Green's functions 
and for the R and L functions allows us to express the 
result for the Green's functions in terms of -Rq+i and L^. 
Hence, we find 

(kll,o;) (A9) 
1 

Now, we are ready to show the expectation value of 
the number current operator vanishes. We can evaluate 
the expectation value of the number current operator in 
Eq. (|18|l by using the Green's functions we have been 
describing above. One finds 

Oq) = a^oQ+i / dw^[Ga+iQ(kll,cj)-Gaa+i(kll,w)] = 0, 
kll 

(AlO) 

which vanishes because the two Green's functions are 
identical in value. Note that this vanishing of the current 
holds for arbitrary size electronic charge reconstruction, 
because it does not involve any linear-response assump- 
tion. Similarly, since the heat-current operator is related 



to the number current operator by a derivative with re- 
spect to time, and that derivative analytically continues 
to an extra power of frequency in the integral over fre- 
quency, we also have that the heat-current operator ex- 
pectation value vanishes, because it involves adding an 
extra power of frequency into the above integrals. This 
then completes the proof that the electronic charge re- 
construction does not have any currents flowing through 
it, so the electric fields that enter the linear-response for- 
malism are the external fields only. 



APPENDIX B: ANALYTIC CONTINUATION OF 
THE FOUR-TIME RESPONSE FUNCTION 
NEEDED FOR THE GENERAL 
JONSON-MAHAN THEOREM PROOF 

We can restrict ourselves to the case n > T2 > T3 > T4 
without loss of generality, because that is the order- 
ing needed to get the relevant correlation functions in 
Eqs. H29|l . First of all, let us introduce a four-time corre- 
lation function (in real time) defined by 

Iabcd (^1,^2,^3,^4) = Iabcd {ti — t,t2 — t,tz — t,t4 — t) 
= {A{ti)B{t2)C{h)D{U)), (Bl) 

where all operators are written in the Heisenberg repre- 
sentation 0{t) = exp[i(H - ^Ji^^)t]0 ey.^[-i{rL - ^iM)t] 
{h = 1). Its Fourier transform gives a four-time spectral 
density (with the constraint o^i + 0^2 + ^^^3 + ^^^4 = due 
to time-translation invariance of equilibrium correlation 
functions) 



+00 -t-00 +00 



lABCD(UJl,UJ2,UJ3,UJ4j = / 2?^ / 2?^ / "^*5Xp[ia;ifi -|-«a;2i2 +«t^3^3 + «t^4t4jjylSCZ3(il,i2,t3,*4) 

— OQ —00 —00 

= ^ XI e"^''AaBifCfpDpid{eu + u}i)S{eif + uj2)S{efp + W3), (B2) 

ilfp 

where is the energy eigenvalue of the quantum many-body state \i) {i. e., the eigenvalue of the operator Ti. — fiJV), 
en satisfies en = £i — On = {i\0\l) is the matrix element of the operator O between the states i and I, and 
Z = e~'^^' is the partition function. The second line in Eq. (jB2p follows from the Lehmann representation by 
inserting appropriate sets of complete states. The spectral density satisfies the following cyclic permutation identities 

lABCD{i^l, (^2,1^3,^^4) = /bcDa(w2, ^3, W4,a;i)e''"i = /(7_DAb(w3, ^^4, ^1, ^^2)6^'^"'+"''^ = /£)ASc(^4, t^l, ^2, t^3, )e"'''^'* 

(B3) 

and transforms under Hermitian conjugation as 

/ABCn(wi,t^2, W3,CJ4) = [/£,tctBtAt(-^^4, -^^3, -t^2, -^l)]*, (B4) 
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with the dagger indicating Hermitian conjugation of the associated operator. Now, the generaUzed function in Eq. H28|) 
can be defined for ti > T2 > > T4 in terms of a generahzed spectral density as 



hoo +00 -|-oo 



-Fa/3(n > T2 > Ta > T4) = J duJi j duj2 J duJsIapiuJl , 0J2, (^3, (^4) exp[-UJlTi - LO2T2 - UJ3T3 - UJ4T4] , (B5) 



-00 —00 



where we introduce the total spectral density in terms of the partial one in Eq. ljB2|l 



zGplaiic jGplaiic 



(B6) 



- I^f c ct c«^, ■('^i''^2,t^3,W4) - /^t ^ _^,.^t (a;i,a;2,W3,W4) ■ 
Next, we can calculate the polarization operators in Eq. (|29|) by taking the appropriate limits and derivatives. Namely, 



+ 00 -\-oo +c>o 



-CSO —00 —00 





for 




= 11 




for 




= 12 


UJl) 


for 


ij 


= 21 


Wl)(w4 - CJ3) 


for 


■ij 


= 22 



(B7) 



Here, one notes that 074 = —lui — lu2 — due to the constraint. Then we replace ivi by + iO^ and construct the 
transport coefficients on the real axis. Finally we take the limit ^ to obtain the dc response 



1 



\-oc +00 



= lini ^[iya/3('^ + iO+) - Lija/3(i^ - iO+)] = 7r/3 / duj2 / dw4/„0(-W2, ^^2, -^^4, ^4)^3 \ (B8) 



-00 —00 



where, by using the identities in Eqs. ljB3|) and (|B4I) . the spectral density in Eq. (|B6p can be reduced to the following 
expression 

Iaf3{-UJ2,UJ2,-UJi,UJ4) ^ -2aH:^^^^tjf3^^ Re /^t^^^^^^^t^^^^^^_(-u;2,W2,-t^4,W4) 

i £ plane j £ plane 

^^^^t {-UJ2,UJ2,-UJi,UJi) . (B9) 



Eqs. HB8|I and (|B9|I are the generalization of the Jonson- 
Mahan theorem to nanostructures; the integrands for the 
charge-charge, heat-charge, charge-heat, and heat-heat 
current operator correlation functions are all related by 
powers of frequency. One can also check that the Onsager 
reciprocal relation holds, where L12 is equal to L2i- This 
follows by using the symmetry relations, and then inter- 
changing the dummy integration variables uj2 and 074. 



and 



+ ^2,-^4,^4) 



= ^c,,,,ci,J^4)/„,,tX^4)e'^-5(c.2-C.4), (Bll) 

respectively, where 



In the case where we neglect vertex corrections, the 
spectral densities in Eq. (|B9|I are equal to 



(-tJ2,t^2, ~'^4,'^4) 



^<=«c^,,(^4)I„,t^^/^4)e^"M(^2 - ^4) 



(BIO) 



■/(cj)ImG'sA(w) 



(B12) 



are the single-particle spectral densities, and we obtain 
the same result as in Eq. H38() . In the general case, when 
vertex corrections are included, one can find the spec- 
tral densities in Eq. I)B2(I from the multitime tempera- 
ture (Matsubara) Green's function in Eq. (|28l) by em- 
ploying spectral relations for the multitime correlation 
functions^^. The full derivation is complex and lengthy. 
In the end it provides no new information with relation 
to the Jonson-Mahan theorem, only an explicit formula 
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for the charge conductivity matrix. Hence, we do not go 
through the derivation here. 
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